Correlated Initial Conditions in Directed Percolation 
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o\ : 

^>«^ , We investigate the influence of correlated initial conditions on the temporal evolution of a 

(d+l)-dimensional critical directed percolation process. Generating initial states with correlations 
{siSi+r) ~ r''~'^ we observe that the density of active sites in Monte-Carlo simulations evolves as 
^ , pit) ~ t". The exponent k depends continuously on a and varies in the range —f3lv\\ < k. < rj. Our 

■ numerical results are confirmed by an exact field-theoretical renormalization group calculation. 

O ' PACS numbers: 05.70.Ln, 64.60.Ak, 64.60.Ht 
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Q . I. INTRODUCTION 

It is well known that initial conditions influence the temporal evolution of nonequilibrium systems. The systems' 
"memory" for the initial state usually depends on the dynamical rules. For example, stochastic processes with a finite 
^ ' temporal correlation length relax to their stationary state in an exponentially short time. An interesting situation 
c/2 ' emerges when a system undergoes a nonequilibrium phase transition where the temporal correlation length diverges. 
This raises the question whether it is possible construct initial states that affect the entire temporal evolution of such 
systems. 

^ To address this question, we consider the example of Directed Percolation (DP) which is the canonical universality 

' class for nonequilibrium phase transitions from an active phase into an absorbing state Q. DP is used as a model 
describing the spreading of some non-conserved agent and may be interpreted as a time-dependent stochastic process 
O , in which particles produce offspring and self-annihilate. Depending on the rates for offspring production and self- 
annihilation such models display a continuous phase transition from a fluctuating active phase into an absorbing state 
without particles from where the system cannot escape. Close to the phase transition the emerging critical behavior is 
characterized by a particle distribution with fractal properties and long-range correlations. The DP phase transition 
is extremely robust with respect to the microscopic details of the dynamical rules |^,^ and takes place even in 1+1 
dimensions. 
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p^j , Monte-Carlo (MC) simulations of critical models with absorbing states usually employ two different types of initial 
' conditions. On the one hand random initial conditions (Poisson distributions) are used to study the relaxation of an 
initial state with a finite particle density towards the absorbing state. In this case the particle density p(t) decreases 
on the infinite lattice asymptotically as (for the definition of the DP scaling exponents P,v±,i'\\,z see Ref. |^) 



Pit) ^ t-^/'^^K (1) 



On the other hand, in so-called dynamic MC simulations ^ , each run starts with a single particle as a localized active 
seed from where a cluster originates. Although many of these clusters survive for only a short time, the number of 
r-| particles n{t) averaged over many independent runs increases as 

8: {nit)}-t+\ (2) 



where rj = — 2j3)/v\\. These two cases seem to represent extremal situations where the average particle number 
either decreases or increases. 

' A crossover between these two extremal cases takes place in a critical DP process that starts from a random initial 
condition at very low density. Here the particles are initially separated by empty intervals of a certain typical size 
wherefore the average particle number first increases according to Eq. (^. Later, when the growing clusters begin to 
interact with each other, the system crosses over to the algebraic decay of Eq. (^ - a phenomenon which is referred 
to as the "critical initial slip" of nonequilibrium systems Q . 

In the present work we investigate whether it is possible to interpolate continuously between the two extremal cases. 
As will be shown, one can in fact generate certain initial states in a way that the particle density on the infinite lattice 
varies as 

p{t) ^ r (3) 
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with a continuously adjustable exponent k in the range 



< K < (4) 

To this end we construct artificial initial configurations with algebraic long-range correlations of the form 

C(r) = (s, s,+r) - r-^-^-^) , (5) 

where () denotes the average over many independent realizations, d the spatial dimension, and Si = 0, f inactive 
and active sites. The exponent cr is a free parameter and can be varied continuously between and 1. The limit of 
long-range correlations cr — > d corresponds to a constant particle density and thus we expect Eq. (|]) to hold. On 
the other hand, the short-range limit ct — » represents an initial state where active sites are separated by infinitely 
large intervals so that the particle density should increase according to Eq. (||). In between we expect p{t) to vary 
algebraically according to Eq. (|^) with an exponent k depending continuously on cr. Our aim is to investigate the 
functional dependence of K{a). 

The effect of power-law correlated initial conditions {<j){0)(t){r)) ~ r^i'^^'^) in case of a quench to the ordered phase 
of systems with nonconserved order parameter has been investigated some time ago by Bray et. al. Such systems 
are characterized by coarsening domains that grow with time as t^^"^ . An important example is the (2-|-l)-dimensional 
Glauber-Ising model quenched to zero temperature. It was observed that long-range correlations are relevant only if 
cr exceeds a critical value ac- Furthermore, it was shown that the relevant regime is characterized by a continuously 
changing exponent in the autocorrelation function A(t) = [(/()(r, <)0(r, 0)] ~ i-('^-<^)/4^ whereas the usual short-range 
scaling exponents could be recovered below the threshold. The results were found to be in agreement with the 
simulation results for the two-dimensional Ising model quenched from T — Tc to T — 0. 

The DP process - the prototype of models with a phase transition from an active phase into an absorbing state - 
is different from coarsening processes. Instead of growing domains the DP process generates fractal clusters of active 
sites with a coherence length ^± which grows as t^^^ where z = v\\/va_- Thus the scaling forms assumed in Ref. |^ 
are no longer valid in the present case. In addition, the field-theoretical description of DP involves nontrivial loop 
corrections and thus we are interested to find out to what extent the results are different from those in Ref. Our 
investigation also sheds some light on the relation between the observed phenomena for correlated initial states, the 
critical initial slip, and scaling laws in time dependent simulations. 

In the present work we focus on the following aspects of the problem: In Sect. || we describe in detail how correlated 
initial states can be constructed in one dimension. Using these states we then perform MC simulations in order to 



numerically estimate the exponent k as a function of cr (see Sect. ID and Fig. ^). In Sect. IV we present a field- 
theoretical renormalization group calculation which generalizes recent results obtained by Wijland et. al. 0. Because 
of a special property of the vertex diagrams and the loop diagrams for the initial particle distribution it is possible to 
derive an exact scaling relation, leading to our main result 

{rj for a < (Tc 

(6) 
^{d- (T ~ (3/v±_) for (T > cr^ 

with the critical threshold cr^ — P/v^. Because of the scaling relation 77 = (i^j^ — 2l3)/v\\ this function is continuous at 
cr = (Jc- The theoretical result is found to be in agreement with our simulation results in one spatial dimension. In 
Sect. we compare the correlations of our constructed initial states with the "natural" correlations that are generated 



by the DP process itself. Finally we summarize our conclusions in Sect. VI. 



II. CONSTRUCTION OF CORRELATED INITIAL STATES 



The construction of artificial correlated particle distributions on a lattice is a highly nontrivial task since the 
lattice spacing and finite size effects lead to deviations that strongly affect the accuracy of the numerical simulations. 
Therefore one has to carefully verify the correlation exponent and fractal dimension of the generated distribution. 
In this Section we describe in detail how such particle distributions can be generated and tested. For simplicity we 
restrict ourselves to initial states in one spatial dimension. 

Let us consider a particle distribution on the real line where particles are separated by empty intervals of length £. 
We assume that these intervals are uncorrelated and power-law distributed according to 

P{C) ^r°' . 1< a < 2 (7) 
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FIG. 2. Numerical verification of the fractal dimension df and the correlation exponent cr on a lattice with 2^^ sites averaged 
over 20000 samples, (a): Number of active boxes m versus the total number of boxes n for various values of a. (b): Estimate 
of the corresponding fractal dimension df. (c)-(d): Analogous estimation of the correlation exponent cr as a function of a. 



This distribution corresponds to a simple fractal set with the fractal dimension df — a — 1, hence the range of a 
is restricted by 1 < a < 2. On a lattice, however, the lattice spacing and the system size L have to be taken into 
account as lower and upper cutoffs for the distribution P{£). The quality of a lattice approximation depends on the 
actual implementation of these cutoffs. It turns out that the accuracy of MC simulations depends strongly on the 
quality of the initial states and therefore the proper implementation of the cutoffs is crucial in the present problem. 

We find that a good approximation is obtained when an (almost) perfect fractal set is projected onto the lattice in 
a way that site i becomes active if at least one element of the fractal belongs to the interval [i,i + 1]. The resulting 
lattice configuration is the minimal set of boxes on the lattice that is needed to cover the fractal set. This projection 
can be efficiently realized on a computer by generating a sequence of points x on the real line separated by intervals 
distributed according to Eq. (|^) with a very small cutoff A^m ^ 1 and projecting it onto the lattice by the following 
prescription (see Fig. |l|): 

1. Start with the empty lattice Si = (i=l,. . . ,L) and let a; be a real variable with the initial value a; = 1. 

2. Let i be the maximal integer for which i < x and turn site i into the active state by setting Si ~ I. 

3. Let Amax = L — x he the current upper cutoff and generate a random number < y < 1 from a flat distribution. 
If y < {Amax / ■Amin)^~°' the construction of the initial state is finished, otherwise continue. 



ill d licit U.10l)liUUl)lUll ill tllC llltCl VcXi 1^^ ^ ^ -/ ^ ^ 

by ^ = z^/''^""-', and continue at step 2. 



4. Generate another random number z from a flat distribution in the interval A^J^ < z < A^^" and increment x 



Notice that step 3 takes the upper cutoff into account by finishing the loop when the generated interval £ would 
exceed the remaining size of the chain L — x. The lower cutoff is processed in step 4 by truncating the allowed range 
of z. 

In order to verify the quality of this approximation, we numerically estimate the fractal dimension dj of the generated 
initial states by box counting. To this end we divide the whole system into n boxes and count the number m of boxes 
that contain at least one active site, averaging over many independent realizations. In Fig. |^a m is plotted against n 
in a double-logarithmic scale. The straight lines indicate that the "true" fractal is well approximated. From the slopes 
we estimate the fractal dimension d/ which is shown in Fig. ^ as a function of a. We also measure the two-point 
correlations in the generated states which should be precisely those of Eq. (||) with cr = a — 1. This can be proven by 
assuming that the intervals are uncorrelated and evaluating a geometric series of the Laplace transform of Pit). In 
order to verify this relation, we estimated <T(a) numerically in Fig. 

In both measurements we find a fairly good agreement with the exact results (dashed lines in Fig. ^). It turns out 
the deviations close to a = 1 can be reduced by increasing the system size while the deviations close to a = 2 are due 
to the lattice spacing and the lower cutoff Amin- 

It should be emphasized that these artificial initial states have a vanishing particle density in the limit L — > oo. 
On a finite lattice, however, a finite density is generated which depends on a and may vary over several decades. By 
increasing the lattice size we therefore reduce the initial particle density which leads to a higher statistical error in 
the subsequent DP process. Thus the optimal system size has to be determined by balancing discretization errors of 
the initial states against statistical errors of the DP process. 



III. NUMERICAL RESULTS 



The time dependent simulations have been performed by using a Domany-Kinzel stochastic cellular automaton 
with parallel updates and periodic boundary conditions. The Domany-Kinzel model is controlled by two parameters 
Pi and p2 and has a whole phase transition line where a DP transition takes place. We performed simulations for 
three different realizations, namely 



• pi = 0.8092, p2 = 0: Wolfram-18 transition, 

• pi = P2 = 0.705485; site percolation, 

• pi = 0.644700, p2 = 0.873762: bond percolation. 
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FIG. 3. Temporal dependence of p for different initial correlations characterized by a = 1...2 in steps of 1/32 (from bottom 
to top). The simulation has been performed on a periodic lattice with L = 8192 sites, a lower cutoff Amin = 0.0001, averaging 
statistically over 1600 samples. 

The results turned out to be the same in all cases within numerical accuracy, although the bond and site percolation 
showed a better scaling law behavior. An efficient, multi-spin coded program has been employed that simulates 32 
replicas for different values a G (1,2), exploiting the advantage that the same random numbers can be used on each 
replica at a given site. This allows us to update 32 systems in parallel stored as an integer vector of length L. The 
lattice size has been varied between L = 128 and 10000, while the number of independent samples was 10000 — 500 
respectively. Finite size effects turned out to be negligible for L > 4096. As shown in Fig. || the quality of the 
power-law for p(t) is quite convincing and extends over the first three decades in the case of a system with L — 8192 
sites. For finite systems all these curves will eventually cross over to the decay after very long time. The slopes 

of the lines in the log-log plot - measuring the exponent K{a) - have been extrapolated from the last two decades 
by standard linear regression analysis. For each a we also determine the corresponding correlation exponent cr(a) 
as explained in the previous Section (see Fig. H). In Fig. ^ we plot K{a) versus a{a). Using the numerically known 
estimates for the DP exponents 

/3 = 0.2765, i/|| = 1.7338, i^_L = 1-0968, z = 1.5807, ry = 0.3137, (8) 

our results are in a fairly good agreement with the theoretical prediction of Eq. (||) (solid line). The deviations for 
small a originate from the lattice cutoff and could be reduced by further increasing the computational effort. 

IV. RENORMALIZATION GROUP CALCULATION 

In this Section we derive Eq. (^) by an exact RG calculation. To this end the DP Langevin equation Q has to be 
extended by an additional term 

dtpix, t) = ap{x, t) - bp^{x, t) -f DV^p{x, t) + r]{x, t) -f Ai0(a;)(5(t) , (9) 

where 4>{x) represents the initial particle distribution and ri{x, t) is a p-dependent noise field with correlations 
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FIG. 4. Estimates for the exponents K(cr) obtained from the slopes of the hnes in Fig. 
by the sohd hne. The dashed hne indicates the "natural" correlations of DP (see Sect. K 



The theoretical prediction is shown 



{ri{x, t)r]{x' , t')) — 7 p{x, t) 5{x — x) S(t — t') . 



(10) 



Although in principle we could directly analyze the Langevin equation, it is often more convenient to derive the 
corresponding effective field theory by introducing a second-quantized bosonic operator representation |]l0| . The 
Langevin equation it then transformed into an effective action S = Sq + Sint + Sipd, where Sq, Smti and Sipd denote 
the free part, the non-linear interaction, and the contribution for the initial particle distribution. Following the 
notation of Ref . H , the respective parts are given by: 



d'^x dt tp{x, t) [dt + A(ct - V^)] i)(x, t) 



<S'mt[V', ^] = y d'^xdtgi}:{x, t) '0(x, t) [i/'(x, t) - ip{x, t)] 

/oo 
d'^xY.^j^'i^^Omx). 
i=i 



(11) 
(12) 
(13) 



The part Sq + Sint is just the usual action of Reggeon field theory [[1 iy2|] . The additional contribution Sipd couples the 
field ^{x,0) with the initial particle distribution (f>{x). It is written in its most general form containing contributions 
of all orders in with independent coefficients Aj. The lowest order contribution in Eq. ( |l3| ) corresponds to 

the term Ai(l>(x)S{t) in the Langevin equation (^. The higher order terms for j > 2 are included because they may 
be generated under RG transformations. However, as we will see below, it is actually sufficient to consider the first 
two contributions. 

The relaxation of a DP process with random initial conditions (pix) = po to its stationary state has been studied 
recently by Wijland et. al. Q using Wilson's dynamic renormalization group approach. At the critical dimension 
dc = 4 the constant contribution Aiipcj) was shown to be relevant while the "Poissonian" contribution A2'!/'^0 turned 
out to be marginal. For d < dc, however, fluctuations corrections cause the coefficient A2 to vanish under RG 
transformations. This allowed the authors to express the so-called critical initial slip exponent 77 [9' in their notation) 
by an exact scaling relation. As mentioned in the Introduction, this exponent describes the initial short-time behavior 
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of p{t) until the correlations generated by the dynamical process are longer than the typical size of empty intervals 
in the initial state |5| so that the system crosses over to the usual decay of Eq. (||). In the present case, however, the 
interval sizes of the initial state are power-law distributed and affect not only the initial but in principle the whole 
temporal evolution. Nevertheless we can use the formalism of Ref. in order to determine the exponent k as a 
function of a. The only difference is that we introduce a field 4>{x) that carries a non-trivial scaling dimension d^. 
Let us first consider a scaling transformation 

x~^x' ^bx, t^t' = bH, (14) 



where z — v±/i'\\ is the anisotropy exponent of DP. Under this transformation the fields in Eqs. ([1 1\)- ( [13D change 
according to 

^ (p'ix') = b-'^'*'(f>{x) 
ip{x, t) il}'{x', t') = b~'^^i}{x, t) (15) 

i!{x, t) -> 4,'{x\ t') = b-'^^i}{x, t) . 

As in Reggeon field theory, the fields V': V' carry the same scaling dimension d^p = d^. The scaling dimension of the 
initial state d^ is related to the fractal dimension df as follows. The number of active sites on a lattice with L 
sites grows as ~ L'^'f . On the other hand, A^ ~ {4')L'^j where (0) denotes the average density of active sites, i.e., 
{(j)) ~ L'^f^'^. Under rescaling this relation turns into {(p)b~'^* = {Lb)'^f~'^, hence 

dcj, = d — df = d — a. (16) 

Under rescaling the contribution Stpd changes by 

„ CXD 

S^pd - Sl^, = y '^'^ E b"-^"^--"^ A, V^Xx) , (17) 
i.e., the coefficients Aj transform according to 



Aj a;. = b'^-^'^^-'^* Aj . (18) 
Defining the anomalous scaling dimensions ry^, 77^ of the fields ip,4>\yy 

di; = ■^{d + rj^) , d^ = ^{d + r]^) (19) 

the scaling dimensions of the coefficients Aj in the mean field approximation is given by 



A, 2 



d{l - j) - - V4> ■ (20) 



This result is expected to hold at the critical dimension dc — 4. Notice that the relevance of the contributions Aj 
decreases with increasing j. In systems with less than four spatial dimensions, fluctuation corrections have to be 
taken into account. Assuming that contributions with j > 3 in Eq. ( [T^ ) are irrelevant, these corrections have been 
computed in Ref. in a one-loop approximation. It was shown that in = 4 — e dimensions the coefficients Ai and 
A2 change under infinitesimal scaling 6 = 1 -I- ^ according to 

rfAj ^ ^^ -d-2^,-,u _ ^2j^ _ ,i^).5£!fl + 0(.') , (22) 

where denotes the surface area of the unit sphere in four dimensions divided by (27r)'* and A is the ultraviolet 
cutoff in momentum space. 

Let us first consider the renormalization of A2. As already noticed in Ref. Q, the diagrams in Eq. ( p^ which are 
linear in A2 are formally identical with those for the renormalization of the nonlinear vertex g to all orders in e (see 
Fig. ^). This observation plays a crucial role in the present problem: Since the renormalization of the vertex g is 
given by 
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FIG. 5. One-loop diagrams for the nonlinear vertex g (left) and the Poissonian contribution A2 (right). The diagrams are 
formally equivalent. The same applies to higher loop diagrams that are linear in A2 and the corresponding vertex diagrams 
(not shown here). Hence the loop corrections of these quantities are are identical to all orders in e. 



dl 



3^ _ 2g^g4 
2 A2 



O(e^) 



(23) 



we may replace the diagrams of the form /S.2[~2g^K/^/\^ + 0(6^)] in Eq. ( p2| ) by the scaling part /S.2[z — d/2 — 3?7^/2] 
of Eq. (G3), resulting in 



dl 



z]+iA2 



-0(6) 



The linear part of this equation is now exact to all orders in e. Since 



z < d^ — z < 



(24) 



(25) 



we arrive at the conclusion that A2 scales to zero for arbitrary initial states in d < 4 spatial dimensions. Thus Eq. ( |2l| ) 
reduces to its scaling part and hence the scaling dimension of Ai is given by 



dAi = = d — d^ 



dd, = cr - d,i 



(26) 



Depending on the sign of , the loop corrections for the initial particle distribution are relevant for a > ac, irrelevant 
for (7 < CTc, and marginal for a = (Tc- 

Continuously changing exponents usually appear when we modify a model with a marginal parameter which is 
invariant under RG transformations. For example, in case of models with infinitely many absorbing states the initial 
density plays that role |^^ . In our case the two-point correlation of the field ■0(a;,t) is characterized by = d — cr 
for r > rc(t) and is left invariant under the coarse-graining RG transformation, where rc(t) — > 00 is a growing spatial 
scale. In the infinite time limit this scaling crosses over to d^ ~ The marginality of the operator Sipd requires 

a — d^ — Q and hence in the infinite time limit we arrive at 



d^{t 



(27) 



In the regime a > ac the exponent n is related to the scaling dimension c?Ai as follows. The particle density is 
expected to vary in the long time limit as 



p{t) c^pot", 



(28) 



where po oc Ai is the initial particle density. Under rescaling ( p^ ) the density p{t) [which is essentially the average 
of tp{x,t)] scales as p{t) h~'''*p{t) while the initial density transforms according to po ^ b^'^^Po- Thus scaling 
invariance of Eq. (|2|) requires that —d^ = zk — d^, i.e., 



2z 



(29) 



In the irrelevant regime cr < CTc, however, the initial density po scales to zero under RG transformations which may be 
interpreted as an initial state with non-interacting active seeds leading to an increase according to Eq. (||). Combining 
these results we arrive at Eq. (0) which is exact to all orders in e. 
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V. "NATURAL" CORRELATIONS IN A DP PROCESS 

It is interesting to compare the artificial initial states of Sect. || with "natural" DP scaling states. Scaling invariance 
predicts that a critical DP process, starting from a fully occupied lattice, evolves towards a state with correlations 

C{r) = {s, s,+r) ^ T~<'l-^ , r « e± (30) 

i.e., the "natural" correlations of the critical DP are given by = d — f3/i'±. Interestingly, the corresponding 

exponent k{(t^^) vanishes. This means that these correlations characterize a situation where the critical DP process 
is almost stationary. From the field-theoretical point of view this is not surprising since k is determined solely by the 
scaling dimension of the initial particle distribution. Nevertheless the result is surprising from the physical point of 
view since the empty intervals in our artificial initial states are uncorrelated whereas such correlations may exist in 
DP scaling states. It seems that these correlations are rather weak so that states with uncorrelated intervals can 
be used as an approximation of DP scaling states. This offers an interesting practical application in DP simulations: 
Instead of starting a critical DP process from random initial conditions and simulating over long transients, one could 
use the artificial states of Sect. ^ with a = as an approximation, followed by a short equilibration period to 
reach the 'true' scaling state of DP. 



VI. CONCLUSIONS 



In the present work we numerically studied the temporal evolution of a (l-l-l)-dimensional critical DP process 
starting from artificially generated correlated initial states. These states approximate a simple fractal set in which 
uncorrelated empty intervals are distributed by P{£) ^ It can be shown that such particle distributions are 

characterized by long-range correlations of the form (siSi+r) ~ r°'~^ where a = a — 1. 

The construction of such correlated states is a technically difficult task since lattice spacing and finite-size effects 
strongly influence the quality of the numerical results, especially close to a « 1 and a « 2. In order to minimize these 
errors, we proposed to project an almost perfect fractal set onto the lattice. Using these states as initial conditions we 
determined the density p{t) in a critical DP process by MC simulations. It varies algebraically with an exponent k that 
depends continuously on a. Therefore correlated initial conditions affect the entire temporal evolution of a critical 
DP process. The numerical estimates for k{(t) are in good agreement with theoretical results from a field-theoretical 
RG calculation (see Fig. |^). Discrepancies for small values of a can be further reduced by increasing the numerical 
effort. 

The RG calculation is valid for arbitrary spatial dimensions d < 4 and predicts a critical threshold CTc = P/i'± where 
K(cr) starts to vary linearly between -fry and — Thus our result in Eq. (^ qualitatively reproduces the scenario 
predicted by Bray et. al. in the context of coarsening processes The exact prediction the critical threshold <Tc 
is possible because of the irrelevance of A2 for arbitrary initial conditions which can be proven by using the formal 
equivalence of the loop expansions in g and A2. It is this property that also allows one to express the critical initial slip 
exponent rj by an exact scaling relation r; = {v^ — 2/3)/j/||. It should be emphasized that some DP models with more 
complicated dynamical rules violate this scaling relation as, for example, the two-species spreading model discussed 
in Ref. [Q. In such models the mentioned equivalence of loop expansions in Fig. |^ is no longer valid. Consequently, 
the corresponding slip exponent 77 takes a different value so that Eq. (H) has to be modified appropriately. 
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